function [psf] = PSF_estimation(k_init,u_img,b_img,R)
% Estiate the PSF with call to PSFest_FINAL
% Author: Liu Yunzhuo, 52543417@qq.com
% Date: 10 January 2022

k_size = sqrt(size(k_init(:), 1));
k = k_init(:);

% 将u_img 处理成 u_hat 并保存供优化算法使用
u_hat = generate_block_circ_u(u_img, k_size);
% 将b_img 处理成 b_hat 并保存供优化算法使用
b_hat = b_img(:);
% 将u_hat' * u_hat 处理成 u_u 并保存供优化算法使用
u_u = u_hat' * u_hat;
% 将u_hat' * b_hat 处理成 u_b 并保存供优化算法使用
u_b = u_hat' * b_hat;


% Perform a deconvolution on the image and the noisy image

% options.Method = 'csd';
% options.MaxIter = 2000;
% psf = minFunc(@PSFest,k, options, u_hat,b_hat,u_u,u_b);

% fminunc
options = optimoptions('fminunc',...
    'SpecifyObjectiveGradient',true,...
    'MaxIterations',1000,...
    'Display','iter');

psf = fminunc(@PSFest,k, options, u_hat,b_hat,u_u,u_b,R);


% fmincon
%     'Algorithm','trust-region-reflective',...
% options = optimoptions('fmincon',...
%     'SpecifyObjectiveGradient',true,...
%     'Display','iter');
% A = [];
% b = [];
% Aeq = [];
% beq = [];
% lb = zeros(size(k));
% ub = 2*ones(size(k));
% noncolon = [];
% psf = fmincon(@PSFest,k, A,b,Aeq,beq,lb,ub,noncolon,options, u_hat,b_hat,u_u,u_b,R);

psf = reshape(psf, [k_size, k_size]);
end
